#Joystick model for tanks according to Bakalis et al. (2017)
 
import openseespy.opensees as op
import vfo.vfo as vfo
import os
import numpy as np
import matplotlib.pyplot as plt
import time
import Units as ut
import scipy as scp


op.wipe();


print('Model generation started...')

direc = './ResultsJM'

if(os.path.exists(direc)=='False'):
    os.mkdir(direc);
    
    
R = 13.9*ut.m #Tank Radius in m
H = 16.5*ut.m #Contained Liquid Height in m
nsp = 8 #Number of beams used to modelling the base plate, Bakalis et al. (2017)
ro = 998*ut.kg_m3 #liquid density in kg/m3 
fillr = 0.95 #Pencentage of tank's height filled with liquid 

m = ro*H*np.pi*R**2*fillr

#Compute impulsive and convective masses, heights, and periods according to Calvi & Nascimbene (2023)
n = np.linspace(0,235,235, dtype='int')
gamma = H*fillr/R
vn = np.zeros(len(n))
for i in range(len(n)):
    vn[i] = np.pi*(2*n[i]+1)/2
    
#Compute modified bessel function of the first kind and its derivative   
I1 = np.zeros(len(n))
I1p = np.zeros(len(n))

for i in range(len(n)):
    I1[i] = scp.special.i1(vn[i]/gamma)
    I1p[i] = scp.special.i0(vn[i]/gamma)- scp.special.i1(vn[i]/gamma)/(vn[i]/gamma)

#Impulsive mass
mi = np.sum(2*m*gamma*(I1/(I1p*vn**3))) 
print(mi/m) 
#Impulsive height
hi = H*np.sum((-1)**n*I1*(vn*(-1)**n-1)/(vn**4*I1p))/np.sum((I1/(I1p*vn**3)))
print(hi)
#Coefficient for Impulsive period, Ci, depends on gamma, page 139 in Calvi and Nascimbene (2023)
Ci = 6.2
tw = 10.61*ut.mm  #wall thickness
rom = 998   #liquid density in kg/m3
Es = 210000000000 #Modulus of elasticity of steel, in Pa
#Impulsive period
Ti = Ci*H*(rom/Es)**0.5/(tw/R)**0.5
#Equivalent impulsive stiffness
ki = 4*np.pi**2*mi/Ti**2
#ki = 100000000000
#print(ki)
ksiI = 0.02 #Damping ratio of impulsive component 
#Define damping constant of impulsive component
ci = ksiI*2*(ki*mi)**0.5

#Convective mass
lam1 = 1.8412  #first root of the first derivative of the Bessel function of the first kind
mc = m*2*np.tanh(lam1*gamma)/(lam1*gamma*(lam1**2-1)) 
#Convctive height
hc = H*(1+(1-np.cosh(lam1*gamma))/(lam1*gamma*np.sinh(lam1*gamma)))
#Fundamental Convective period
Tc = 2*np.pi*(R/ut.g)**0.5/(lam1*np.tanh(lam1*gamma))**0.5
#Equivalent convective stiffness
kc = 4*np.pi**2*mc/Tc**2

ksiC = 0.005 #Damping ratio of convective component 
#Define damping constant of convective component
cc = ksiC*2*(kc*mc)**0.5

op.model('basic', '-ndm', 3, '-ndf', 6) 

#Define central node
nodeTagCent = 1
op.node(nodeTagCent,0,0,0)

#Define impulsive mass nodes
nodeTagi1 = 2
nodeTagi2 = 3
op.node(nodeTagi1, 0,0,hi)
op.node(nodeTagi2, 0,0,hi)

#op.equalDOF(nodeTagi1,nodeTagi2,3,4,5,6)

#Define impulsive mass
op.mass(nodeTagi2, mi,mi,0)

#Define convective mass nodes
nodeTagc1 = 4
nodeTagc2 = 5
op.node(nodeTagc1, 0,0,hc)
op.node(nodeTagc2, 0,0,hc)

#op.equalDOF(nodeTagc1,nodeTagc2,3,4,5,6)

#Define convective mass
op.mass(nodeTagc2, mc,mc,0)

#Define nodes on tank circumference
dth = 2*np.pi/nsp
nodeTagCirc1 = 10
nodeTagCirc2 = 100

for i in range(0,nsp):
    op.node(nodeTagCirc1+i, round(R*np.cos(dth*i),5),round(R*np.sin(dth*i),5),0)
    op.node(nodeTagCirc2+i, round(R*np.cos(dth*i),5),round(R*np.sin(dth*i),5),0)
    op.fix(nodeTagCirc2+i, 1,1,1,1,1,1)
    #op.equalDOF(nodeTagCirc2+i,nodeTagCirc1+i,1,2,4,5,6)

#Define Materials

#Parameters of Flag-shape material model used to model the base plate uplift behaviour, calibrated using baseplate model 
matTag1 = 1

ePf1 = 350
ePd1 = 0.0038
k1 = ePf1/ePd1

ePf2 = 550
ePd2 = 0.09
k2 = (ePf2-ePf1)/(ePd2-ePd1)

sigAct = ePf1
beta =0.38
epsSlip = 0.0
epsBear = 0.2

ePf3 = 1800
ePd3 = 0.4

k3 = (ePf3-ePf2)/(ePd3-ePd2)

rBear = k3/k1

#Define material of springs modelling the base plate uplift 
op.uniaxialMaterial('SelfCentering',matTag1,k1,k2,sigAct,beta,epsSlip,epsBear,rBear)

matTag2 = 2
K = 800/0.3
op.uniaxialMaterial('Elastic',matTag2,K,0,ut.Ubig)

#Define impulsive and convective elastic materials
ImpMatTag = 3
op.uniaxialMaterial('Elastic',ImpMatTag,ki)
ConvMatTag = 4
op.uniaxialMaterial('Elastic',ConvMatTag,kc)

#Define impulsive and convective viscous materials
ImpMatTagD = 5
op.uniaxialMaterial('Viscous',ImpMatTagD,ci,1)
ConMatTagD = 6
op.uniaxialMaterial('Viscous',ConMatTagD,cc,1)

rigidMatTag = 7
op.uniaxialMaterial('Elastic',rigidMatTag,ut.Ubig)

flexMatTag = 8
op.uniaxialMaterial('Elastic',flexMatTag,ut.Usmall)


ParMatTag = 9
op.uniaxialMaterial('Parallel',ParMatTag,matTag1,matTag2)

#Define elements
Transf1 = 1
Transf2 = 2
op.geomTransf('Linear', Transf1,1,0,0)  # Geometric transformation for rigid elastic elements 
op.geomTransf('Linear', Transf2,0,0,1) 
#Define rigid elastic elements
Ae = 1000000000
Ee =  210
Ize = 8333333333300000
Iye = 8333333333300000
Ge = Ee/(2*(1+0.3))
Je = 8333.333333300000

eleTagiBeam = 1
eleTagcBeam = 2

op.element('elasticBeamColumn', eleTagiBeam, nodeTagCent, nodeTagi1, Ae, Ee, Ge, Je, Iye, Ize, Transf1)
op.element('elasticBeamColumn', eleTagcBeam, nodeTagCent, nodeTagc1, Ae, Ee, Ge, Je, Iye, Ize, Transf1)

#op.rigidLink('beam',nodeTagi1,nodeTagCent)
#op.rigidLink('beam',nodeTagc1,nodeTagCent)

eleTagi = 3
eleTagc = 4

op.element('zeroLength', eleTagi, nodeTagi1,nodeTagi2,'-mat',ImpMatTag,ImpMatTag,rigidMatTag,rigidMatTag,rigidMatTag,rigidMatTag,'-dir',1,2,3,4,5,6)
op.element('zeroLength', eleTagc, nodeTagc1,nodeTagc2,'-mat',ConvMatTag,ConvMatTag,rigidMatTag,rigidMatTag,rigidMatTag,rigidMatTag,'-dir',1,2,3,4,5,6)

eleTagSpoke = 100
eleTagBP = 1000
for i in range(0,nsp):
    op.element('elasticBeamColumn',eleTagSpoke+i,nodeTagCent,nodeTagCirc1+i, Ae, Ee, Ge, Je, Iye, Ize, Transf2)
    #op.rigidLink('beam',nodeTagCent,nodeTagCirc1+i)
    op.element('zeroLength',eleTagBP+i,nodeTagCirc2+i,nodeTagCirc1+i,'-mat',rigidMatTag,rigidMatTag,ParMatTag,flexMatTag,flexMatTag,flexMatTag,'-dir',1,2,3,4,5,6)


#Run Gravity Analysis

constTS = 1
op.timeSeries('Constant', constTS)
patternTagG = 1
op.pattern('Plain',patternTagG,constTS)
#op.load(nodeTagi2,0,0,0,0,0,0)
#op.load(nodeTagc2,0,0,0,0,0,0)

Tol = 1e-7 # convergence tolerance for test
NstepGravity = 10
DGravity = 1/NstepGravity
op.integrator('LoadControl', DGravity) # determine the next time step for an analysis
op.numberer('RCM') # renumber dof's to minimize band-width (optimization), if you want to
op.system('BandGen') # how to store and solve the system of equations in the analysis
op.constraints('Plain') # how it handles boundary conditions
op.test('NormDispIncr', Tol, 100) # determine if convergence has been achieved at the end of an iteration step
op.algorithm('Newton') # use Newton's solution algorithm: updates tangent stiffness at every iteration
op.analysis('Static') # define type of analysis static or transient
op.analyze(NstepGravity) # apply gravity
op.loadConst('-time', 0.0) #maintain constant gravity loads and reset time to zero
print('Static analysis completed')


omega = []
freq =  []
T = []

lamb = op.eigen('-fullGenLapack',4)

for lam in lamb:
    omega.append((lam)**0.5)
    freq.append((lam)**0.5/(2*np.pi))
    T.append((2*np.pi)/(lam)**0.5)

for i in  range(len(T)):
    print('T'+str(i+1)+' = '+str(T[i])+' s')


#Define Recorder
op.recorder('Node', '-file', direc+'/dispI.out','-time', '-node', nodeTagi2, '-dof', 1, 'disp')
op.recorder('Node', '-file', direc+'/UpliftI.out','-time', '-node', nodeTagCirc1+4,'-dof',3,'disp')
op.recorder('Element', '-file', direc+'/force1.out','-time','-ele', eleTagBP, 'force')
op.recorder('Element', '-file', direc+'/force2.out','-time','-ele', eleTagBP+1, 'force')
op.recorder('Element', '-file', direc+'/force3.out','-time','-ele', eleTagBP+2, 'force')
op.recorder('Element', '-file', direc+'/force4.out','-time','-ele', eleTagBP+3, 'force')
op.recorder('Element', '-file', direc+'/force5.out','-time','-ele', eleTagBP+4, 'force')
op.recorder('Element', '-file', direc+'/force6.out','-time','-ele', eleTagBP+5, 'force')
op.recorder('Element', '-file', direc+'/force7.out','-time','-ele', eleTagBP+6, 'force')
op.recorder('Element', '-file', direc+'/force8.out','-time','-ele', eleTagBP+7, 'force')
op.recorder('Element', '-file', direc+'/forceI.out','-time','-ele', eleTagi, 'force')
op.recorder('Element', '-file', direc+'/forceC.out','-time','-ele', eleTagc, 'force')

vfo.createODB(model="Tank", loadcase="Pushover")

#Run Cyclic Pushover Analysis
print('Running Cyclic Pushover...')
#displist = [2*ut.mm,0,-2*ut.mm,0,4*ut.mm,0,-4*ut.mm,0,6*ut.mm,0,-6*ut.mm,0,9*ut.mm,0,-9*ut.mm,0,12*ut.mm,0,-12*ut.mm,0,23*ut.mm,0,-23*ut.mm,0,35*ut.mm,0,-35*ut.mm,0,46*ut.mm,0,-46*ut.mm,0,58*ut.mm,0,-58*ut.mm,0,70*ut.mm,0,-70*ut.mm,0,90*ut.mm,0,-90*ut.mm,0,105*ut.mm,0,-105*ut.mm,0,115*ut.mm,0,-115*ut.mm,0,115*ut.mm,0,-115*ut.mm,0]
displist = [100*ut.mm]
dstep = 0.01*ut.mm
IDctrlNode = nodeTagi2          # node where disp is read for disp control
IDctrlDOF  = 1                  # degree of freedom read for disp control (1 = x displacement)

testParams = [1.e-4, 50]

op.numberer('RCM') # renumber dof's to minimize band-width (optimization), if you want to
op.system('FullGeneral') # how to store and solve the system of equations in the analysis
op.constraints('Lagrange') # how it handles boundary conditions
op.test('EnergyIncr',*testParams) # determine if convergence has been achieved at the end of an iteration step
op.algorithm('Newton') # use Newton's solution algorithm: updates tangent stiffness at every iteration
op.analysis('Static') # define type of analysis static or transient

POtag = 100;
linTS = 2
op.timeSeries('Linear',linTS)

op.pattern('Plain', POtag, linTS)     
op.load(IDctrlNode,1,0,0,0,0,0)
op.reactions()
force=[op.nodeReaction(nodeTagCirc2,1)]
ok = 0
for j in range(0, len(displist)):
    if ok != 0:
        break
    dispnew = displist[j]
    if j > 0:
        dispold = displist[j-1]
    else:
        dispold = 0
    disp = dispnew - dispold
    if disp > 0:
        dsteps = dstep
    else:
        dsteps = -dstep
    nstep = int(disp/dsteps)
    op.integrator('DisplacementControl', IDctrlNode, IDctrlDOF, dsteps) # determine the next time step for an analysis
    for i in range(nstep):
        ok = op.analyze(1)
        
        if ok != 0:
            testParams2 = [1.e-4, 2000]
            op.test('EnergyIncr',*testParams2)
            op.algorithm('Newton','-initial')
            print("Trying Newton with Initial Tangent ..")
            ok = op.analyze(1)
            op.test('NormDispIncr',*testParams) 
            op.algorithm('Newton') 
        
        if ok != 0:
            op.algorithm('Broyden',50)
            op.test('RelativeNormDispIncr',*testParams2)
            print("Trying Broyden ..")
            ok = op.analyze(1) 
            op.test('NormDispIncr',*testParams) 
            op.algorithm('Newton') 
            
        if ok != 0:
            op.algorithm('NewtonLineSearch')
            print("Trying NewtonWithLineSearch ..")
            ok = op.analyze(1)
            op.algorithm('Newton') # use Newton's solution algorithm: updates tangent stiffness at every iteration
            
        if ok == 0:
            force.append(op.nodeReaction(nodeTagi1,1))
            dsp = round(op.nodeDisp(IDctrlNode,IDctrlDOF),3)
            #print(f'Displacement {dsp} reached of {dispnew}')
        
        else:
            break

if(ok != 0):
    print('Problem with Pushover')
else:
    print('Done')
#print(force)        
#ani = vfo.animate_deformedshape(model="Tank", loadcase="Pushover")
vfo.plot_model(show_nodes = "yes",show_nodetags="yes")
#vfo.plot_modeshape(modenumber=1,scale=100)
#Define parameters of Pinching4 model for the base springs
